<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/1999/REC-html401-19991224/loose.dtd">
<html lang="en">
<head>
	<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
	<title>Simulating B-mode Ultrasound Images Example (k-Wave)</title>
	<link rel="stylesheet" href="kwavehelpstyle.css" type="text/css">
	<meta name="description" content="Simulating B-mode Ultrasound Images Example.">
</head>

<body><div class="content">

<h1>Simulating B-mode Ultrasound Images Example</h1>

<p>This example illustrates how k-Wave can be used for the simulation of B-mode ultrasound images (including tissue harmonic imaging) analogous to those produced by a modern diagnostic ultrasound scanner. It builds on the <a href="example_us_defining_transducer.html">Defining An Ultrasound Transducer</a>, <a href="example_us_beam_patterns.html">Simulating Ultrasound Beam Patterns</a>, and  <a href="example_us_transducer_as_sensor.html">Using An Ultrasound Transducer As A Sensor</a> examples.</p>

<p>
    <ul>
        <li><a href="matlab:edit([getkWavePath('examples') 'example_us_bmode_linear_transducer.m']);" target="_top">Open the file in the MATLAB Editor</a></li>
        <li><a href="matlab:run([getkWavePath('examples') 'example_us_bmode_linear_transducer']);" target="_top">Run the file in MATLAB</a></li>
    </ul>
</p>

<p>Note, this example generates a B-mode ultrasound image from a 3D scattering phantom using <code><a href="kspaceFirstOrder3D.html">kspaceFirstOrder3D</a></code>. Compared to ray-tracing or Field II, this approach is very general. In particular, it accounts for nonlinearity, multiple scattering, power law acoustic absorption, and a finite beam width in the elevation direction. However, it is also computationally intensive. Using a modern GPU and the Parallel Computing Toolbox (with <code>'DataCast'</code> set to <code>'gpuArray-single'</code>),  each scan line takes around 3 minutes to compute. Using a modern desktop CPU (with <code>'DataCast'</code> set to <code>'single'</code>), this increases to around 30 minutes. In this example, the final image is constructed using 96 scan lines. This makes the total computational time around 4.5 hours using a single GPU, or 2 days using a CPU.</p>

<p>To allow the simulated scan line data to be processed multiple times with different settings, the simulated RF data is saved to disk. This can be reloaded by setting <code>run_simulation = false</code> within the example m-file. The data can also be downloaded from <a href="http://www.k-wave.org/datasets/example_us_bmode_scan_lines.mat">http://www.k-wave.org/datasets/example_us_bmode_scan_lines.mat</a></p>

<h2>Contents</h2>
<div>
	<ul>
		<li><a href="#heading2">Simulation approach</a></li>
		<li><a href="#heading3">Defining the scattering medium</a></li>
		<li><a href="#heading4">Running the simulation</a></li>
		<li><a href="#heading5">Processing steps</a></li>
		<li><a href="#heading6">Fundamental and harmonic imaging</a></li>
		<li><a href="#heading7">Extending the simulations</a></li>		
	</ul>
</div>

<a name="heading2"></a>
<h2>Simulation approach</h2>

<p>To minimise the computational overhead of running multiple 3D simulations with large grid sizes, the total number of elements in the ultrasound transducer is set to the number of active elements (in this case 32). This allows a smaller computational grid to be used. For each new scan line, the values for <code>medium.sound_speed</code> and <code>medium.density</code> within the computational grid are updated based on a larger map of the scattering phantom (see schematic below). This has the effect of simulating a larger transducer in which the active elements are swept across the transducer.</p>

<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_01.png" style="width:648px;height:414px;" alt="">

<p>Note, the same result could be achieved by allocating the complete transducer and then progressively updating <code>transducer.active_elements</code>. However, this requires a larger computational grid for the simulation of each scan line.</p>

<a name="heading3"></a>
<h2>Defining the scattering medium</h2>

<p>The contrast in an ultrasound image is generated by acoustic reflections from changes in impedance (sound speed and density). Within the human body, these heterogeneities occur across multiple length scales. At a macroscopic level, different tissue types such as fat or muscle have different bulk material properties. Specular and diffractive scattering from heterogeneities on this scale manifest as noticeable edges in the ultrasound image. At a microscopic level, there is also a difference in the acoustic properties between the cells and the extra cellular matrix in which they preside. Diffusive scattering from heterogeneities on this scale ultimately results in the speckle pattern that is characteristic of ultrasound images. (The speckle pattern is created by the constructive and destructive interference of spatial and temporal variations in these reflections. The exact speckle pattern is also strongly dependent on the characteristics of the ultrasound probe.)</p>

<p>To account for these different effects, the ultrasound phantom must similarly be heterogeneous across multiple scales. In this example, a scattering phantom is created by modulating the mean sound speed and density at each grid point using random Gaussian noise. To simulate bulk tissue contrast, three spherical regions with increased scattering and impedance are defined using <code><a href="makeBall.html">makeBall</a></code>. The centers of the three spheres are co-aligned with center of the transducer in the elevation direction. A schematic of the simulation layout and a slice through the scattering phantom in the beam plane are shown below.</p>

<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_02.png" style="width:560px;height:420px;" alt="">
<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_03.png" style="width:560px;height:420px;" alt="">

<a name="heading4"></a>
<h2>Running the simulation</h2>

<p>The simulation is run in a similar way to the previous ultrasound examples. The defined transducer is used to replace both the <code>source</code> and <code>sensor</code> inputs. The scan lines are then simulated sequentially. For each scan line, the values for <code>medium.sound_speed</code> and <code>medium.density</code> within the computational domain are updated based on the larger scattering phantom. The pressure at each active transducer element (in this case 32) is then returned and stored.</p>

<pre class="codeinput">
<span class="comment">% loop through the scan lines</span>
for scan_line_index = 1:number_scan_lines

    <span class="comment">% load the current section of the medium</span>
    medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :);
    medium.density = density_map(:, medium_position:medium_position + Ny - 1, :);
        
    <span class="comment">% run the simulation</span>
    sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});

    <span class="comment">% extract the scan line from the sensor data</span>
    scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
        
    <span class="comment">% update medium position</span>
    medium_position = medium_position + transducer.element_width;

end
</pre>

<a name="heading5"></a>
<h2>Processing steps</h2>

<p>After simulation, a number of steps are required to process the raw transducer data into an ultrasound image. In this example, the processing follows the basic steps performed by a modern diagnostic ultrasound scanner. These are illustrated in the block diagram below. More detailed information on each of these steps can be found in any reference book on ultrasound imaging, for example, "Diagnostic Ultrasound: Principles and Instruments" by Frederick Kremkau, or "Diagnostic Ultrasound Imaging" by Thomas Szabo.</p>

<p></p>
<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_04.png" style="width:673px;height:103px;" alt="">
<p></p>

<p><b>1. Receive Beamforming</b></p>
<p>First, the received pressure signals at each active transducer element are processed into a single scan line. This is done using the <code>scan_line</code> method of the <code><a href="kWaveTransducer.html">kWaveTransducer</a></code> class. The beamforming delays and weights are automatically calculated based on the transducer focus distance and apodization settings.</p>

<pre class="codeinput">
<span class="comment">% extract the scan line from the sensor data</span>
scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
</pre>

<p>After beamforming, the first part of each scan line is also windowed to remove the input signal from the recorded sensor data. This step is required because k-Wave allows the transducer elements to both send and receive at the same time. (This is not physically possible using a standard ultrasound probe.)</p>

<p></p>
<p><b>2. Time Gain Compensation</b></p>
<p>The acoustic waves generated by the ultrasound transducer lose energy as they propagate through the tissue due to acoustic attenuation (absorption and scattering) and focussing effects. Consequently, reflections from deeper tissue features will appear weaker in the recorded signals. To compensate for this, a simple correction is applied assuming the attenuation within tissue follows the relation:</p>
<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_eq_01.png" style="width:107px;height:23px;" alt="">
<p>where alpha is the attenuation coefficient in dB/(MHz cm), f is the frequency in MHz (chosen here as the transmit frequency), and d is the total (round-trip) distance travelled in cm. Here the distance variable is created assuming that t = 0 corresponds to the middle of the transmitted tone burst.</p>

<pre class="codeinput">
<span class="comment">% create radius variable assuming that t0 corresponds to the middle of the
% input signal</span>
t0 = length(input_signal) * kgrid.dt / 2;
r = c0 * ( (1:length(kgrid.t_array)) * kgrid.dt / 2 - t0);    <span class="comment">% [m]</span>

<span class="comment">% define absorption value and convert to correct units</span>
tgc_alpha_db_cm = medium.alpha_coeff * (tone_burst_freq * 1e-6)^medium.alpha_power;
tgc_alpha_np_m = tgc_alpha_db_cm / 8.686 * 100;

<span class="comment">% create time gain compensation function based on attenuation value and
% round trip distance</span>
tgc = exp(tgc_alpha_np_m * 2 * r);

<span class="comment">% apply the time gain compensation to each of the scan lines</span>
scan_lines = bsxfun(@times, tgc, scan_lines);
</pre>

<p></p>
<p><b>3. Frequency Filtering</b></p>
<p>The scan lines are then filtered using a Gaussian filter centered about the transmit frequency or the second harmonic. This reduces the effects of noise outside the transmit frequency range, or in the latter case, obtains the harmonic signal for tissue harmonic imaging. The filtering is performed using <code><a href="gaussianFilter.html">gaussianFilter</a></code> by specifying the required filter center frequency and percentage bandwidth.</p>

<pre class="codeinput">
<span class="comment">% filter the scan lines using both the transmit frequency and the second harmonic</span>
scan_lines_fund = gaussianFilter(scan_lines, 1/kgrid.dt, tone_burst_freq, 100, true);
scan_lines_harm = gaussianFilter(scan_lines, 1/kgrid.dt, 2 * tone_burst_freq, 30, true);
</pre>

<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_05.png" style="width:560px;height:420px;" alt="">
<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_06.png" style="width:560px;height:420px;" alt="">

<p></p>
<p><b>4. Envelop Detection</b></p>
<p>The high frequency variations in the scan lines are then removed in a demodulation process known as envelope detection. Here, this is achieved using <code><a href="envelopeDetection.html">envelopeDetection</a></code> which extracts the signal envelope using the Hilbert transform.</p>

<pre class="codeinput">
<span class="comment">% envelope detection</span>
scan_lines_fund = envelopeDetection(scan_lines_fund);
scan_lines_harm = envelopeDetection(scan_lines_harm);
</pre>

<p></p>
<p><b>5. Log Compression</b></p>
<p>The dynamic range of the processed scan lines often exceeds the gray scale range that can be displayed by a computer monitor or visually perceived. To reduce this dynamic range, the scan lines are log compressed using <code><a href="logCompression.html">logCompression</a></code> with a normalised compression ratio of 3.</p>

<pre class="codeinput">
<span class="comment">% normalised log compression</span>
compression_ratio = 3;
scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
scan_lines_harm = logCompression(scan_lines_harm, compression_ratio, true);
</pre>

<p></p>
<p><b>6. Scan Conversion</b></p>
<p>Finally, the processed scan lines are spatially remapped to give an image resolution suitable for display. In this example, the generated image is upsampled by a factor of 2 using <code><a href="matlab: doc interp2">interp2</a></code>.</p>

<pre class="codeinput">
<span class="comment">% upsample the image using linear interpolation</span>
scale_factor = 2;
scan_lines_fund = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
scan_lines_harm = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
</pre>

<p>A plot of the central scan line after each processing step is given below.</p>

<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_07.png" style="width:560px;height:420px;" alt="">

<a name="heading6"></a>
<h2>Fundamental and harmonic imaging</h2>

<p>The resulting B-mode images are shown below. The images look realistic and contain many of the features and artifacts seen in diagnostic images from commercial ultrasound scanners, for example, speckle and shadowing. The harmonic image has improved resolution and contrast.</p>

<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_08.png" style="width:560px;height:420px;" alt="">
<img vspace="5" hspace="5" src="images/example_us_bmode_linear_transducer_09.png" style="width:560px;height:420px;" alt="">

<a name="heading7"></a>
<h2>Extending the simulations</h2>

<p>This example can be used as a template for the simulation of a wide range of B-mode ultrasound images. For instance, the scattering phantom could be replaced with one generated using 3D MRI or CT data. It is also possible to simulate the effect of different signal processing techniques, for example, frequency compounding (where the images produced using different band pass filters are averaged), spatial compounding (where the images produced when the transducer is slightly rotated are averaged), and superharmonic imaging (where the higher frequency harmonics are used to form the image).</p>

</div></body></html>